{
 "cells": [
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:12:22.528287Z",
     "start_time": "2025-07-31T04:12:22.116682Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.linalg import lstsq\n",
    "path = \"D:/7_30data/\""
   ],
   "id": "8b8ec8ff1986deb3",
   "outputs": [],
   "execution_count": 1
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:12:22.828576Z",
     "start_time": "2025-07-31T04:12:22.813470Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def iterative_refinement_solve(V, L_rho, max_iter=15, tolerance=1e-15, verbose=True):\n",
    "    \"\"\"\n",
    "    迭代精化求解复数线性系统 S @ a = b，其中 S = V^H @ V, b = V^H @ L_rho\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    V : ndarray, shape (m, n)\n",
    "        输入复数矩阵\n",
    "    L_rho : ndarray, shape (m,)\n",
    "        右端复数向量\n",
    "    max_iter : int\n",
    "        最大迭代次数\n",
    "    tolerance : float\n",
    "        收敛容差\n",
    "    verbose : bool\n",
    "        是否输出详细信息\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    a : ndarray\n",
    "        解向量 (复数形式)\n",
    "    \"\"\"\n",
    "\n",
    "    # 构建正规方程 S @ a = b\n",
    "    S = V.conj().T @ V  # V^H @ V\n",
    "    b = V.conj().T @ L_rho  # V^H @ L_rho\n",
    "\n",
    "    # 计算并输出矩阵的秩信息\n",
    "    rank_S = np.linalg.matrix_rank(S)\n",
    "    augmented = np.hstack((S, b.reshape(-1, 1)))\n",
    "    rank_augmented = np.linalg.matrix_rank(augmented)\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"rank(S) = {rank_S}\")\n",
    "        print(f\"rank([S | b]) = {rank_augmented}\")\n",
    "\n",
    "        # 输出条件数信息\n",
    "        cond_S = np.linalg.cond(S)\n",
    "        print(f\"Condition number of S: {cond_S:.4e}\")\n",
    "        print(\"\\n=== Iterative Refinement Process ===\")\n",
    "\n",
    "    # 检查系统的相容性\n",
    "    if rank_S != rank_augmented:\n",
    "        print(\"Warning: System is inconsistent (rank(S) ≠ rank([S|b]))\")\n",
    "        print(\"Using least squares solution...\")\n",
    "\n",
    "    # 步骤1: 获得初始解 (使用最小二乘法)\n",
    "    try:\n",
    "        a = lstsq(S, b, lapack_driver='gelsy')[0]\n",
    "        if verbose:\n",
    "            print(\"Initial solution computed using least squares\")\n",
    "    except Exception as e:\n",
    "        print(f\"Failed to compute initial solution: {e}\")\n",
    "        return np.zeros(S.shape[1], dtype=complex)\n",
    "\n",
    "    # 初始残差\n",
    "    initial_residual = np.linalg.norm(S @ a - b)\n",
    "    if verbose:\n",
    "        print(f\"Initial residual norm: {initial_residual:.4e}\")\n",
    "\n",
    "    # 迭代精化过程 - 记录残差历史用于更智能的停止判断\n",
    "    residual_history = []\n",
    "\n",
    "    for i in range(max_iter):\n",
    "        # 步骤2: 计算当前残差 r = b - S @ a\n",
    "        r = b - S @ a\n",
    "        residual_norm = np.linalg.norm(r)\n",
    "        residual_history.append(residual_norm)\n",
    "\n",
    "        if verbose:\n",
    "            print(f\"Iteration {i+1}: Residual norm = {residual_norm:.4e}\")\n",
    "\n",
    "        # 步骤3: 检查收敛\n",
    "        if residual_norm < tolerance:\n",
    "            if verbose:\n",
    "                print(f\"Converged after {i+1} iterations\")\n",
    "            break\n",
    "\n",
    "        # 步骤4: 求解修正方程 S @ da = r\n",
    "        try:\n",
    "            da = lstsq(S, r, lapack_driver='gelsy')[0]\n",
    "        except Exception as e:\n",
    "            if verbose:\n",
    "                print(f\"Iteration {i+1} failed: {e}\")\n",
    "            break\n",
    "\n",
    "        # 步骤5: 更新解 a = a + da\n",
    "        a = a + da\n",
    "\n",
    "        # 更智能的停滞检测：只有在残差连续增加或完全停滞时才停止\n",
    "        if i >= 5:  # 至少迭代5次再考虑停止\n",
    "            recent_residuals = residual_history[-4:]  # 最近4次的残差\n",
    "\n",
    "            # 检查是否连续3次残差都没有显著改善\n",
    "            if len(recent_residuals) >= 4:\n",
    "                improvements = []\n",
    "                for j in range(1, len(recent_residuals)):\n",
    "                    if recent_residuals[j-1] > 0:\n",
    "                        improvement = (recent_residuals[j-1] - recent_residuals[j]) / recent_residuals[j-1]\n",
    "                        improvements.append(improvement)\n",
    "\n",
    "                # 如果连续3次改善都小于1%，且残差仍然远大于机器精度，则可能真的停滞了\n",
    "                if (len(improvements) >= 3 and\n",
    "                    all(imp < 0.01 for imp in improvements[-3:]) and\n",
    "                    residual_norm > 1e-14):  # 只有远大于机器精度时才考虑停滞\n",
    "\n",
    "                    if verbose:\n",
    "                        print(f\"Convergence stalled after {i+1} iterations (minimal improvement)\")\n",
    "                    break\n",
    "\n",
    "    # 计算最终残差\n",
    "    final_residual = np.linalg.norm(S @ a - b)\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"Final residual norm for S @ a - b: {final_residual:.4e}\")\n",
    "\n",
    "        # 输出改进情况\n",
    "        improvement_factor = initial_residual / final_residual if final_residual > 0 else float('inf')\n",
    "        print(f\"Improvement factor: {improvement_factor:.2e}\")\n",
    "\n",
    "        # 检查解的性质\n",
    "        max_imag = np.max(np.abs(np.imag(a)))\n",
    "        if max_imag > 1e-12:\n",
    "            print(f\"Warning: Solution has significant imaginary parts: {max_imag:.4e}\")\n",
    "\n",
    "    # 输出必要的信息\n",
    "    print(f\"rank(S) = {rank_S}\")\n",
    "    print(f\"rank([S | b]) = {rank_augmented}\")\n",
    "    print(f\"Residual norm for S @ a - b: {final_residual:.4e}\")\n",
    "\n",
    "    return a"
   ],
   "id": "2dcf31c2eb74bac1",
   "outputs": [],
   "execution_count": 2
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:12:22.911755Z",
     "start_time": "2025-07-31T04:12:22.867411Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def iterative_refinement_solve_real_output(V, L_rho, max_iter=15, tolerance=1e-15, verbose=True):\n",
    "    \"\"\"\n",
    "    迭代精化求解，返回实数解 (如果虚部可忽略)\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    V : ndarray, shape (m, n)\n",
    "        输入复数矩阵\n",
    "    L_rho : ndarray, shape (m,)\n",
    "        右端复数向量\n",
    "    max_iter : int\n",
    "        最大迭代次数\n",
    "    tolerance : float\n",
    "        收敛容差\n",
    "    verbose : bool\n",
    "        是否输出详细信息\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    a : ndarray\n",
    "        解向量 (如果虚部可忽略则返回实数，否则返回复数)\n",
    "    \"\"\"\n",
    "\n",
    "    # 调用复数版本\n",
    "    a_complex = iterative_refinement_solve(V, L_rho, max_iter, tolerance, verbose)\n",
    "\n",
    "    # 检查是否可以转换为实数\n",
    "    max_imag = np.max(np.abs(np.imag(a_complex)))\n",
    "\n",
    "    if max_imag < 1e-12:  # 虚部可忽略\n",
    "        if verbose:\n",
    "            print(\"Converting to real solution (imaginary parts negligible)\")\n",
    "        return np.real(a_complex)\n",
    "    else:\n",
    "        if verbose:\n",
    "            print(f\"Keeping complex solution (max imaginary part: {max_imag:.4e})\")\n",
    "        return a_complex\n",
    "\n",
    "# 高精度版本 - 完全按照您原始成功的算法\n",
    "def high_precision_iterative_solve(V, L_rho, max_iter=20, tolerance=1e-16, verbose=True):\n",
    "    \"\"\"\n",
    "    高精度迭代精化求解 - 完全复制您原始成功的算法逻辑\n",
    "    不使用过早的停滞检测，让算法充分迭代到机器精度\n",
    "    \"\"\"\n",
    "\n",
    "    # 构建正规方程 S @ a = b\n",
    "    S = V.conj().T @ V  # V^H @ V\n",
    "    b = V.conj().T @ L_rho  # V^H @ L_rho\n",
    "\n",
    "    # 计算并输出矩阵的秩信息\n",
    "    rank_S = np.linalg.matrix_rank(S)\n",
    "    augmented = np.hstack((S, b.reshape(-1, 1)))\n",
    "    rank_augmented = np.linalg.matrix_rank(augmented)\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"rank(S) = {rank_S}\")\n",
    "        print(f\"rank([S | b]) = {rank_augmented}\")\n",
    "\n",
    "    # 步骤1: 获得初始解 (使用最小二乘法)\n",
    "    a = lstsq(S, b, lapack_driver='gelsy')[0]\n",
    "    if verbose:\n",
    "        print(\"Initial solution computed\")\n",
    "\n",
    "    # 迭代精化过程 - 简单版本，不做复杂的停滞检测\n",
    "    for i in range(max_iter):\n",
    "        # 计算当前残差\n",
    "        r = b - S @ a\n",
    "        residual_norm = np.linalg.norm(r)\n",
    "\n",
    "        if verbose:\n",
    "            print(f\"Iteration {i+1}: Residual norm = {residual_norm:.4e}\")\n",
    "\n",
    "        # 简单的收敛检查\n",
    "        if residual_norm < tolerance:\n",
    "            if verbose:\n",
    "                print(f\"Converged to tolerance after {i+1} iterations\")\n",
    "            break\n",
    "\n",
    "        # 求解修正方程并更新\n",
    "        try:\n",
    "            da = lstsq(S, r, lapack_driver='gelsy')[0]\n",
    "            a = a + da\n",
    "        except:\n",
    "            if verbose:\n",
    "                print(\"Iterative refinement failed, using current solution\")\n",
    "            break\n",
    "\n",
    "    # 计算最终残差\n",
    "    final_residual = np.linalg.norm(S @ a - b)\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"Final residual norm for S @ a - b: {final_residual:.4e}\")\n",
    "\n",
    "    # 输出必要的信息 (与您的需求格式一致)\n",
    "    print(f\"rank(S) = {rank_S}\")\n",
    "    print(f\"rank([S | b]) = {rank_augmented}\")\n",
    "    print(f\"Residual norm for S @ a - b: {final_residual:.4e}\")\n",
    "\n",
    "    return a"
   ],
   "id": "8284d209691253ac",
   "outputs": [],
   "execution_count": 3
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:12:22.973764Z",
     "start_time": "2025-07-31T04:12:22.968247Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 推荐使用的函数 - 直接替换您原来的函数\n",
    "def improved_iterative_solve(V, L_rho):\n",
    "    \"\"\"\n",
    "    推荐版本：使用高精度迭代精化求解\n",
    "    完全复制您原始成功算法的逻辑，应该能达到 2e-16 的精度\n",
    "    \"\"\"\n",
    "    solution = high_precision_iterative_solve(V, L_rho, max_iter=20, tolerance=1e-16, verbose=True)\n",
    "\n",
    "    # 检查是否可以转换为实数\n",
    "    max_imag = np.max(np.abs(np.imag(solution)))\n",
    "    if max_imag < 1e-12:\n",
    "        return np.real(solution)\n",
    "    else:\n",
    "        return solution"
   ],
   "id": "481fd5b9d4359029",
   "outputs": [],
   "execution_count": 4
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:21:17.516380Z",
     "start_time": "2025-07-31T04:21:14.633745Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 使用示例_1\n",
    "V = np.load(path+'V.npy').astype(np.complex128)\n",
    "L_rho = np.load(path+'L_rho.npy').astype(np.complex128)\n",
    "\n",
    "#调用函数\n",
    "solution = improved_iterative_solve(V, L_rho)\n",
    "\n",
    " #或者使用详细版本\n",
    "solution = iterative_refinement_solve(V, L_rho, max_iter=20, tolerance=1e-16)\n",
    "\n",
    "pass"
   ],
   "id": "83ee4887e03306b6",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Initial solution computed\n",
      "Iteration 1: Residual norm = 1.4462e-10\n",
      "Iteration 2: Residual norm = 2.0580e-11\n",
      "Iteration 3: Residual norm = 8.2318e-11\n",
      "Iteration 4: Residual norm = 2.0580e-11\n",
      "Iteration 5: Residual norm = 2.0580e-11\n",
      "Iteration 6: Residual norm = 2.0580e-11\n",
      "Iteration 7: Residual norm = 2.1084e-16\n",
      "Iteration 8: Residual norm = 2.1328e-16\n",
      "Iteration 9: Residual norm = 2.0580e-11\n",
      "Iteration 10: Residual norm = 2.0580e-11\n",
      "Iteration 11: Residual norm = 1.9571e-16\n",
      "Iteration 12: Residual norm = 2.0580e-11\n",
      "Iteration 13: Residual norm = 6.1739e-11\n",
      "Iteration 14: Residual norm = 2.0580e-11\n",
      "Iteration 15: Residual norm = 2.0580e-11\n",
      "Iteration 16: Residual norm = 2.0580e-11\n",
      "Iteration 17: Residual norm = 2.0580e-11\n",
      "Iteration 18: Residual norm = 2.0580e-11\n",
      "Iteration 19: Residual norm = 2.0553e-16\n",
      "Iteration 20: Residual norm = 2.0580e-11\n",
      "Final residual norm for S @ a - b: 4.1159e-11\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Residual norm for S @ a - b: 4.1159e-11\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Condition number of S: 9.7122e+80\n",
      "\n",
      "=== Iterative Refinement Process ===\n",
      "Initial solution computed using least squares\n",
      "Initial residual norm: 1.4462e-10\n",
      "Iteration 1: Residual norm = 1.4462e-10\n",
      "Iteration 2: Residual norm = 2.0580e-11\n",
      "Iteration 3: Residual norm = 8.2318e-11\n",
      "Iteration 4: Residual norm = 2.0580e-11\n",
      "Iteration 5: Residual norm = 2.0580e-11\n",
      "Iteration 6: Residual norm = 2.0580e-11\n",
      "Iteration 7: Residual norm = 2.1084e-16\n",
      "Iteration 8: Residual norm = 2.1328e-16\n",
      "Iteration 9: Residual norm = 2.0580e-11\n",
      "Iteration 10: Residual norm = 2.0580e-11\n",
      "Convergence stalled after 10 iterations (minimal improvement)\n",
      "Final residual norm for S @ a - b: 1.9571e-16\n",
      "Improvement factor: 7.39e+05\n",
      "rank(S) = 2\n",
      "rank([S | b]) = 2\n",
      "Residual norm for S @ a - b: 1.9571e-16\n"
     ]
    }
   ],
   "execution_count": 9
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:24:03.013813Z",
     "start_time": "2025-07-31T04:24:01.307687Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 使用示例_2\n",
    "np.random.seed(42)\n",
    "V = np.random.randn(4, 273) + 1j * np.random.randn(4, 273)\n",
    "L_rho = np.random.randn(4) + 1j * np.random.randn(4)\n",
    "solution = iterative_refinement_solve(V, L_rho, max_iter=20, tolerance=1e-16)\n",
    "\n",
    "# 使用示例_3\n",
    "np.random.seed(100)\n",
    "V = np.random.randn(4, 273) + 1j * np.random.randn(4, 273)\n",
    "L_rho = np.random.randn(4) + 1j * np.random.randn(4)\n",
    "solution = iterative_refinement_solve(V, L_rho, max_iter=20, tolerance=1e-16)"
   ],
   "id": "36038a8a90d88acc",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Condition number of S: 1.8539e+19\n",
      "\n",
      "=== Iterative Refinement Process ===\n",
      "Initial solution computed using least squares\n",
      "Initial residual norm: 2.8866e-14\n",
      "Iteration 1: Residual norm = 2.8866e-14\n",
      "Iteration 2: Residual norm = 1.7601e-14\n",
      "Iteration 3: Residual norm = 1.7712e-14\n",
      "Iteration 4: Residual norm = 1.7279e-14\n",
      "Iteration 5: Residual norm = 1.8575e-14\n",
      "Iteration 6: Residual norm = 1.9667e-14\n",
      "Iteration 7: Residual norm = 1.8971e-14\n",
      "Iteration 8: Residual norm = 1.8019e-14\n",
      "Iteration 9: Residual norm = 1.7507e-14\n",
      "Iteration 10: Residual norm = 1.8729e-14\n",
      "Iteration 11: Residual norm = 1.8044e-14\n",
      "Iteration 12: Residual norm = 1.7860e-14\n",
      "Iteration 13: Residual norm = 1.7866e-14\n",
      "Iteration 14: Residual norm = 1.8976e-14\n",
      "Iteration 15: Residual norm = 1.8901e-14\n",
      "Convergence stalled after 15 iterations (minimal improvement)\n",
      "Final residual norm for S @ a - b: 1.7803e-14\n",
      "Improvement factor: 1.62e+00\n",
      "Warning: Solution has significant imaginary parts: 4.5904e-02\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.7803e-14\n"
     ]
    }
   ],
   "execution_count": 10
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-07-31T04:27:02.934170Z",
     "start_time": "2025-07-31T04:27:01.476196Z"
    }
   },
   "cell_type": "code",
   "source": "",
   "id": "de6e0a50d4050357",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Condition number of S: 1.0462e+19\n",
      "\n",
      "=== Iterative Refinement Process ===\n",
      "Initial solution computed using least squares\n",
      "Initial residual norm: 4.5320e-14\n",
      "Iteration 1: Residual norm = 4.5320e-14\n",
      "Iteration 2: Residual norm = 1.4481e-14\n",
      "Iteration 3: Residual norm = 1.6181e-14\n",
      "Iteration 4: Residual norm = 1.5660e-14\n",
      "Iteration 5: Residual norm = 1.6769e-14\n",
      "Iteration 6: Residual norm = 1.5028e-14\n",
      "Iteration 7: Residual norm = 1.5818e-14\n",
      "Iteration 8: Residual norm = 1.5589e-14\n",
      "Iteration 9: Residual norm = 1.5697e-14\n",
      "Iteration 10: Residual norm = 1.5809e-14\n",
      "Iteration 11: Residual norm = 1.5670e-14\n",
      "Convergence stalled after 11 iterations (minimal improvement)\n",
      "Final residual norm for S @ a - b: 1.7517e-14\n",
      "Improvement factor: 2.59e+00\n",
      "Warning: Solution has significant imaginary parts: 3.3175e-02\n",
      "rank(S) = 4\n",
      "rank([S | b]) = 4\n",
      "Residual norm for S @ a - b: 1.7517e-14\n"
     ]
    }
   ],
   "execution_count": 11
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}
